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We calculate analytically the critical connectivity A" c of Random Threshold Networks (RTN) 
for homogeneous and inhomogeneous thresholds, and confirm the results by numerical simulations. 
We find a super-linear increase of K c with the (average) absolute threshold \h\, which approaches 
if c (|/i|) ~ /i 2 /(21n \h\) for large \h\, and show that this asymptotic scaling is universal for RTN with 
Poissonian distributed connectivity and threshold distributions with a variance that grows slower 
than h 2 . Interestingly, we find that inhomogeneous distribution of thresholds leads to increased 
propagation of perturbations for sparsely connected networks, while for densely connected networks 
damage is reduced; the cross-over point yields a novel, characteristic connectivity Kd, that has no 
counterpart in Boolean networks. Last, local correlations between node thresholds and in-degree 
are introduced. Here, numerical simulations show that even weak (anti-)correlations can lead to 
a transition from ordered to chaotic dynamics, and vice versa. It is shown that the naive mean- 
field assumption typical for the annealed approximation leads to false predictions in this case, since 
correlations between thresholds and out-degree that emerge as a side-effect strongly modify damage 
propagation behavior. 



PACS numbers: 05.45.-a, 05.65.+b, 89.75.-k, 89. 75. Da 



I. INTRODUCTION 



Many systems in nature, technology and society can be 
described as complex networks with some flow of matter, 
energy or information between the entities the system is 
composed of; examples are neural networks, gene regu- 
latory networks, food webs, power grids and friendship 
networks. Often, in particular when the networks con- 
sidered are very large, many details of the topological 
structure as well as of the dynamical interactions be- 
tween units are unknown, hence, statistical methods have 
to be applied to gain insight into the global properties 
of these systems. In this spirit, Kauffman [l[ Q intro- 
duced the notion of Random Boolean Networks (RBN), 
originally as a simplified model of gene regulatory net- 
works (GRN). In a RBN of size N, each node i receives 
inputs from < k < N other nodes (with k usually 
either considered to be constant, or distributed accord- 
ing to a Poissonian with average K <C N), and updates 
its state according to a Boolean function /j of its in- 
puts; the subscript i indicates that Boolean functions 
vary from site to site, usually assigned at random to 
each node. It was shown that RBN exhibit a percolation 
transition from ordered to chaotic dynamics at a critical 
connectivity K = K c = 2. Since interactions in RBN 
are asymmetric and hence a Hamiltonian does not exist, 
mean-field techniques have to be applied for analytical 
calculation of critical points, for example the so-called 
annealed approximation (annealed approximation) intro- 
duced by Derrida and Pomeau [1, H[ Hj] . In the annealed 
approximation, random perturbations are applied to ini- 
tial dynamical states, and random ensemble techniques 
are applied to determine whether the so- induced "dam- 
age" spreads over the network or not. Recent research 
has revealed many surprising details of RBN dynamics at 
criticality, e.g. super-polynomial scaling of the number 



of different dynamical attractors (fixed points or peri- 
odic cycles) with N [7| (while Kauffman assumed it to 
scale ~ y/~N as well as analytically derived scaling 
laws for mean attractor periods [8| and for the number 
of frozen and relevant nodes in RBN 

[1G3. Similarly, it 
was shown recently that dynamics in finite RBN exhibits 
considerable deviations from the annealed a ppr oximation 
(that is exact only in the limit N — > oo) pHll2|. Boolean 
network models have been applied successfully to model 
the dynamics of real biological systems, e.g. the seg- 
ment polarity network of Drosophila (l3| . dynamics and 
robustness of the yeast cell cycle network [14]], damage 
spreading in knock-out experiments [151 ] as well as estab- 
lishment of position information [TfJ ] and cell differentia- 
tion [l7j in development. Other models explicitly evolve 
RBN topology according to local rewiring rules coupled 
to local order parameters of network dynamics (e.g., the 
local rate of state changes) , and investigate the resulting 
self-organized critical state [H, [Til H3| ■ 



A drawback of RBN is the fact that, in spite of their 
discrete nature (which makes them easy to simulate on 
the computer in principle), the time needed to compute 
their dynamics in many instances scales exponential in N 
and K, and often large statistical ensembles are needed 
for unbiased statistics due to the strongly non-ergodic 
character [2l[ of RBN dynamics. For this reason, there 
exists considerable interest in simplified models of RBN 
dynamics, as, for example, Random Threshold Networks 
(RTN), that constitute a subset of RBN. 

In RTN, states of network nodes are updated accord- 
ing to a weighted sum of their inputs plus a threshold h, 
while interaction weights take (often discrete and binary) 
positive or negative values assigned at random. The crit- 
ical connectivity, calculated by means of the annealed 
approximation, was found to deviate slightly from RBN 
[23 . 123I HH ; this analysis was extended to RTN dynamics 
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including stochastic update errors (25j. In particular, it 
was found that phase transitions in RTN with scale- free 
topologies [111, [26| substantially differ from both RTN 
with homogeneous or Poissonian distributed connectiv- 
ity and scale- free RBN [27j . Further, dynamics in finite 
RTN with k = const. = 2 inputs per node recently was 
found to be surprisingly ordered, including, e.g., globally 
synchronized oscillations (28|. Other approaches, that 
apply learning algorithms as well as ensemble techniques, 
present evidence that information processing of static [29[ 
or time-variant (30j external inputs is optimized at criti- 
cality in both RBN and RTN. 

In this paper, we extend the theoretical analysis of 
RTN in a number of respects. First, we calculate the crit- 
ical connectivity K c for arbitrary thresholds h < 0, and 
generalize this derivation for the first time to inhomoge- 
neously distributed thresholds hi that can vary from node 
to node. This generalization, that introduces an addi- 
tional level of complexity to RTN dynamics, is motivated 
by recent observations of strong variations in regulatory 
dynamics from gene to gene in real GRN, caused by, for 
example, the frequent occurrence of canalizing functions 
(2l| and the abundance of regulatory RNA in multicel- 
lular organisms which strongly influence the expression 
levels and -patterns of (regulatory) proteins [3l|. Using 
the annealed approximation and additional approxima- 
tion techniques, we derive a general scaling relationship 
between critical connectivity K c and (average) absolute 
node threshold \h\, and show that K c (\h\) asymptotically 
approaches a unique scaling law K c (\h\) ~ h 2 /(2 In \h\) 
for large \h\. Evidence is presented that this asymp- 
totic scaling law is universal for RTN with Poissonian 
distributed connectivity and threshold distributions with 
a variance that grows slower than \h\ 2 . Convergence 
against this scaling law is rather slow (logarithmic in \h\); 
we show that, for finite \h\, scaling behavior can be ap- 
proximated well locally by power laws if c (|/i|) ~ \h\ a 
with 3/2 < a < 2. 

Further, we establish that damage propagation func- 
tions of RTN with homogeneous thresholds \h\ and of 
RTN with inhomogeneous thresholds with the same av- 
erage \h\ = \h\ intersect at characteristic connectivities 
Kd(\h\) > K c (\h\), which implies that for K < Kd, ran- 
dom distribution of thresholds tends to increase damage, 
while for K > Kd, the opposite holds. Evidence is pre- 
sented that Kd{\h\) converges to an asymptotic scaling 
law Kd{\h\) ~ h 2 . We compare the scaling of Kd to the 
corresponding case of random Boolean networks (RBN) 
with inhomogeneously distributed bias, parameterized in 
terms of a bias parameter 1/2 < p < 1. It is shown that 
Kd is not defined for RBN in the limit p — * 1, which 
corresponds to — > oo in RTN. Hence, Kd constitutes 
a truly novel, not previously known concept, yielding a 
new characteristic connectivity which is well-defined only 
for RTN. 

Last, we investigate the effect of correlations between 
thresholds hi and in-degree fcj, while keeping all other 
network parameters constant. We find that even small 



positive correlations can induce a transition from super- 
critical (chaotic) to subcritical (ordered) dynamics, while 
anti-correlations have the opposite effect. It is shown 
that the naive mean-field assumption typical for the an- 
nealed approximation leads to false predictions in this 
case. Even in the simplest case, where only in-degree 
and (absolute) threshold are correlated, complete infor- 
mation about topology, including the output side, has to 
enter statistics, and the order of averages becomes im- 
portant. 



II. RANDOM THRESHOLD NETWORKS 

A Random Threshold Network (RTN) consists of N 
randomly interconnected binary sites (spins) with states 
oi = ±1. For each site i, its state at time t+1 is a function 
of the inputs it receives from other spins at time t: 



CTi(*+l) = sgn(/i(t)) 



with 



/»(*) = y^cjj(Tj(t) + ^ 



(i) 



(2) 



where are the interaction weights. If i does not 



receive signals from j, one has Cj 



0, otherwise, 



interaction weights take discrete values c%j — ±1, +1 or 
— 1 with equal probability. In the following discussion we 
assume that the threshold parameter takes integer values 
hi < Further, we define sgn(0) = -1. H The 

N network sites are updated synchronously. Notice that 
we depart from the well-studied case hi — const. = in 
two respects: hi can take arbitrary values hi < 0, and it 
can differ from node to node (inhomogeneous thresholds). 

Let us now have a closer look on network topology. Let 
K be the average connectivity, i.e. the average number 
of inputs (outputs) per site, and let us assume that each 
interaction weight has equal probability p = K/N to 
take a non-zero value. Further, let us consider the limit 
of sparsely connected networks with K <C N. Under 
these assumptions, the statistical distribution pk of in- 
and out-degrees follows a Poissonian: 



Pk 



K k 



-K 



(3) 



Further, we study the case where in- and out-degree 
distributions differ: while the out-degree is still dis- 
tributed according to a Poissonian, the in-degree distri- 
bution exhibits a power-law tail, i.e. 



Pk in oc k" 



(4) 



with 2 < 7 < 4. 
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III. CALCULATING THE CRITICAL LINE 
A. Uniform threshold h < 

We start with the simplest case and assume that all 
network sites have identical integer threshold values hi = 
h < 0. The case h > is not studied here, as it may lead 
to the pathological outcome of nodes set to an active 
state <Ji = +1, though they receive only inhibitory inputs 
Cij < 0. 

Let us first calculate the probability for damage 
spreading p s (k), i.e. the probability that a node with 
k inputs changes its state, if one of its input states is 
flipped. A straight-forward extension of the combinato- 
rial analysis carried out in [2J] for the special case h = 
yields 



P. 



(k,\h\) = fc" 1 • 2-<* +1 > • 



( k 

(k + \h\ + l)- ( k+jhj+i 

V 2 



+(k-\h\ + l) • fe _|M 



l+i 



2 -C*-i) ((*-!) 

z 1 k+\h\-l 

\ 2 



(5) 
(6) 



for odd k — \h\ with k > \h\, and 



P.(k,\h\) 



jfe-i.a-Cfc+i) 

+ (k+\h\ + 2) ■ 
L \ k±\h\ 

V 9 



(k~\h\) 



k-\h\ 



k 

k+\h\+2 



(7) 
(8) 



for even k — \h\ with k > \h\ (for a detailed derivation, 
please refer to appendix A). Notice that Eqs. ([6|) and 
([5]) are similar, yet not identical to the corresponding 
relations derived in [25| for RTN with probabilistic time 
evolution; in particular, for the RTN with deterministic 
dynamics as studied here, the relation p° dd ( k) = p s (k— 1) 
holds only for the special case \h\ =0, whereas for \h\ > 0, 
p s (k) exhibits an oscillatory behavior (Fig. [1]). 

If we know the statistical distribution function pk of 
the in-degree, the average damage spreading probability 
then simply follows as [24| 



N 



( Ps ) = PkPs(k + l,\h\), 
k=\h\ 



(9) 



where (.) indicates the average over the ensemble of all 
possible network topologies that can be generated accord- 
ing to the degree-distribution pk ■ In the case of a Poisson 
distributed connectivity with average degree degree K, it 
follows 



JY 



(p s )(K,\h\) = 



-,-K 



E ^- P s(k + i,\h\). 

k=\h\ 



(10) 
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FIG. 1: Probability p s (k, \h\) of damage propagation, for dif- 
ferent values of the threshold \h\, as a function of the number 
of inputs k. For large k, the curves asymptotically approach 
p s ~ 1/vfe (dashed line). Notice the oscillatory behavior for 
\h\>0. 
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FIG. 2: Expectation value d of damage one time step after a 
one-bit perturbation, as a function of the average connectivity 
K, and different (homogeneous) thresholds \h\ {\h\ = (+), 
\h\ = 1 (X), |A| = 2 (*), \h\ = 3 (□), |A| = 4 (0). Solid 
curves are the corresponding analytical results obtained from 
the annealed approximation. 



Let us now apply the so-called annealed approximation 
0, which averages the effect of perturbations over the 
whole ensemble of possible network topologies and all 
possible state configurations; in this approximation, the 
expected damage d after one update time step, given a 
one-bit perturbation at time t — 1 then follows as 



d(t+l) = (p s )(K,\h\)-K, 



(11) 



where . denotes the average over all possible network 
topologies and all possible state configurations. If we 
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FIG. 3: Average damage d(K) one time step after a one-bit 
perturbation, for Poisson-distributed connectivity with aver- 
age degree K, and Poisson-distributed negative thresholds 
with average absolute value \h\; points are data from numer- 
ical simulations of RTN (ensemble averages over 100000 dif- 
ferent network realizations for each data point), lined curves 
are analytical solutions (annealed approximation) . Numerical 
data where sampled for \h\ = (+), \h\ = 0.3 (X), \h\ = 1.0 
(*), \h\ = 1.5 (squares), \h\ = 2.5 (o), \h\ — 3.5 (triangle) and 
|ft| = 5.0 (+). 

apply a sufficiently large (but finite) upper limit N to 
the sum in Eq. (|10|) . we can numerically evaluate this 
formula with any desired accuracy. Figure shows the 
results for the first five values of negative h of RTN with 
Poissonian distributed connectivity, compared to mea- 
surements obtained from numerical simulations of large 
ensembles of randomly generated instances of RTN, indi- 
cating an excellent match between theory and simulation. 



B. Poisson distributed thresholds 

Let us now consider the more general case of non- 
uniform thresholds, i.e., networks where each site i has 
assigned an individual threshold hi < 0. In the sim- 
plest case, we can imagine that the final thresholds re- 
sulted from iterated, random decrementations (starting 
from h = for all sites), until a certain average threshold 
h is reached - this process results in Poisson distributed 
thresholds hi. If threshold assignment is independent 
from the (also Poisson distributed) in-degree, the prob- 
abilities for k and h simply multiply, and the resulting 
average damage propagation probability is 

_ _ N N i> k \h\\h\ 

( Ps m\h\) = e-^ J2 E ^mrp^+m), 

|ft|=ofc=|ft| ■' 

(12) 

where \h\ is the average absolute threshold. 

Figure [3] demonstrates that the expected damage 



FIG. 4: Comparison of damage spreading in networks with 
homogenenous thresholds \h\ = const, (solid lines, threshold 
values \h\ as indicated) vs. networks with inhomogeneous 
thresholds distributed according to a Poissonian with the 
same average threshold \h\ (curves with data points, \h\ = 1 
(+), \h\ = 2 (x), \h\ = 3 (*) and \h\ = 4 (□); results obtained 
from the annealed approximation. 

dt+i(K, \h\) resulting from a one-bit perturbation at time 
i, as predicted from this annealed approximation over 
both degree- and threshold distribution, exhibits excel- 
lent agreement with the results obtained from numeri- 
cal simulations of randomly generated RTN ensembles. 
It is an interesting question how the dynamics of RTN 
with inhomogeneous thresholds compares to RTN with 
homogeneous thresholds. Figure 2] shows d(K) for RTN 
with different homogeneous \h\ = const, and the corre- 
sponding inhomogeneous RTN with_ Poisson-distributed 
thresholds with the same average \h\ = \h\, as obtained 
from the annealed approximation. One observes that for 
small K, the curves for RTN with inhomogeneously dis- 
tributed thresholds are systematically above those of the 
corresponding homogeneous RTN, i.e., the randomiza- 
tion of node thresholds increases dynamical disorder - 
also, the critical connectivities i^ c (|/i|) (intersections with 
the line d = 1) are shifted to smaller values. However, 
one also realizes that the curves intersect in the super- 
critical phase at characteristic connectivities i.e., 
for K > Kd(\h\), inhomogeneity in thresholds actually 
reduces damage. 

C. Universal scaling of the critical line 

If we again assume a one-bit perturbation at time t, the 
critical line K c (\h\), that separates the ordered and the 
chaotic phase of RTN dynamics, is given by the condition 

d(t + 1) = (p s )(K c (\h\), \h\) ■ K c (\h\) = 1. (13) 

Again, we can apply Eq. (JTUJ) to solve this equation for 
arbitrary h < 0, however, numerical evaluation is almost 
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FIG. 5: Logarithm of the average damage, In [d(K)], as calcu- 
lated from the annealed approximation, for different values of 
\h\ (\h\ = 10 (+), \h\ = 20 (X), \h\ = 40 (*) and \h\ = 60 (□)). 
The corresponding solid curves are obtained from Eq. (I14[l . 
For not to small K, one finds that Eq. (|14fl approximates the 
true damage function very well. 
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FIG. 6: Scaling behavior of the critical connectivity A c (|ft|) as 
a function of the (homogeneous) node threshold \h\, log- log- 
plot. Data points + are solutions obtained from the annealed 
approximation of Eq. Q5) , the solid curve is obtained from 
setting Eq. (I14|l to zero. The dashed line shows the asymp- 
totic scaling behavior stated in Eq. (). 



impossible for \h\ > 80 due to exponentially diverging 
computing time caused by evaluation of the sum in Eq. 
(fT5)l for large A [36| . For estimation of the scaling behav- 
ior of A c (|ft|) for larger |ft|, we are interested in a good 
approximation that does not require summation over the 
whole network topology, and hence neglect the variation 
in k, considering damage propagation in the mean field 
limit k = const, w A (for details, see Appendix B). Us- 
ing the Stirling approximation: n\ « n"e~ ,l v / 27m, this 
leads to the following approximation for the logarithm of 



the damage: 



ln[d(A,|ft|)] 



- { In A- A -In 



|ft|ln 



K+\h\ 



A 



C 



K-\h\_ 

with C = In [y/2/ , K S j ; solving this equation for 
ln[d(A c (|ft|),|ft|] = 



(14) 



(15) 



then yields the critical connectivity A c (|/i|). Figure [S] 
shows that this approximation is very accurate even for 
considerably small, finite \h\. In particular, one can show 
that for | ft | > 10 the relative error e between the approx- 
imation of Eq. (|15p and the result obtained from the 
annealed approximation vanishes ~ |ft| _1 (Fig. [7]). 

Still, Eq. (fT5")) has to be solved numerically to calculate 
A c (|ft|), and hence does not yield information about the 
scaling behavior in the limit \h\ — > oo. A first insight into 
the expected scaling can be obtained from an analysis of 
the scaling behavior of the maximum of p s (k, \h\) with 
respect to |ft|; if we restrict our analysis to even k — |ft|, 
kmax is given by the condition 



Ap s =p s (k,\h\)-p s (k-2,\h\),^0 



(16) 



or, more accurately, we have to find the minimum of the 
absolute value |Ap s /Afc| of the 'discrete derivative' of 
p s (k, \h\) for even k— |ft|, with Afc = const. = 2. Inserting 
Eq. © then yields 



&Ps 



-fe+3 



(fc-3)! 



[(fc+|ft|-3)/2]![(fe-|ft|-3)/2]! 
(k- l)(fc-2) 



<17) 



(k + \h\ + l)(k- \h\ -1) 



Obviously, the pre-factor on the right hand-side is always 
positive; consequently, in order to determine the maxi- 
mum of p s (k, |ft|), we have to solve the equation 



(fc-l)(fc-2) 



(fc+|fc| + l)(*-|fe|-l. 
Using simple algebra, one can show that 

kmax | ft] ~T" 1 



-1 = 0. 



(18) 



(19) 



solves this equation, i.e. the maximum of p s {k, \h\) scales 
quadratically with |ft|. Since p s (k, |ft|) for |ft| 3> van- 
ishes both for small and large k, it is plausible that the 
scaling behavior of A c is dominated by the leading be- 
havior of the maximum of the distribution, i.e. should 
scale ~ /(|ft|)|ft| 2 , where contributions from the tails of 
the distribution are considered in /(|ft|). 

A more detailed analysis carried out in appendix [Cl 
takes into account that, for large A and |ft|, according 
to the central limit theorem the Binomial distribution 
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FIG. 7: Crosses (X): Relative error si between the approxi- 
mation of Eq. (|14[) and the result obtained from the annealed 
approximation, as a function of \h\. For \h\ > 15, ei vanishes 
oc |/t| _1 ; straight line with slope —1 shown for comparison. 
Data points (+): Relative error £2 between the approxima- 
tion of Eq. (|14[) and the asymptotic scaling of Eq. (|21|l ; £2 
goes to zero logarithmically (compare to dashed curve). 
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FIG. 8: Optimal exponents a of power- laws K c « a\h\ a that 
approximate the scaling function K c (\h\), as shown in Fig. [6] 
as a function of \h\. The dashed curves are the corresponding 
asymptotic estimates of Eq. (|23|l and Eq. (I24|l . 



that characterizes the damage propagation function Eq. 
§E§ can be replaced by a Gaussian, consequently, the ex- 
pected damage is approximated very well by 



d(K,\h\) = K 



1 



2vrA' 



exp 



2K 



(20) 



Taking logarithms and inserting into Eq. (|15[) then yields 
the asymptotic scaling 



(21) 



of the critical connectivity K c . Figure [6] demonstrates 
the convergence of the critical line (straight-lined curve 
and data points) against this asymptote (dashed curve). 

For finite \h\, we notice that there are substantial con- 
tributions from additional terms that vanish only loga- 
rithmically, and hence an approximation based on Eq. 
([21) would substantially underestimate K c . This can be 
appreciated clearly from Fig. which demonstrates the 
slow (logarithmic) convergence of the error £2(|/i|) made 
by application of Eq. (|2"Tj) for finite \h\. 

From Fig. [51 it is also evident that, for finite \h\, Eq. 
(j2"Tj) overestimates the slope dK c /d\h\. One can show 
that, for finite \h\, K c (\h\) is better approximated locally 
by power-laws of the form 



K c (\h\) « a(\h\) ■ \hfW 



(22) 



with 3/2 < a < 2. We confirmed this intuition by nu- 
merically inserting candidate solutions with fixed a into 
Eq. HU), and solving for the values of \h\ and a where 
the deviation from the true curve iT c (|/i|) becomes mini- 
mal; inverting this relation, we obtain the optimal power 
law exponents a(|/i|) as a function of \h\ (Fig. [8l for de- 
tails, see appendix E). Again, we can apply the Gaussian 
approximation for the damage propagation function to 
derive upper (lower) bounds for the finite size scaling of 
a(|/i|) and a(|/i|), which yield (cf. appendix E) 



and 



a(\h\)*2- 



a(\h\) 



1 



bx\h\ 



2hx\h\ 



(23) 



(24) 



lim KJ\h\) = —. — j-rr 
|h|— 00 Vl u 2\n\h\ 



Figure [5] shows that the true optimal values are system- 
atically below (a) or above (a) these curves, demonstrat- 
ing the non-trivial scaling behavior of the critical line for 
finite \h\, which is significantly different from the sim- 
ple asymptotic behavior in the thermodynamic limit (Eq. 

Let us now investigate the scaling behavior of K c for 
networks with inhomogeneous thresholds. Figure[H]shows 
that, for finite \h\, the critical line K c (\h\) for RTN 
with inhomogeneous thresholds is always below the corre- 
sponding values for homogeneous J h \ ; the absolute differ- 
ence AK c (\h\ := \K^(\h\) - Ki(\h\ = \h\)\ between both 
curves, however, increases only linearly in with \h\ (inset 
of Fig. [21), where K^(\h\) is the critical connectivity for 
homogeneous \h\, and A*(|ft,|) the corresponding value for 
inhomogeneously distributed \h\ with mean \h\ = \h\. 

Intuitively, this is straight-forward to understand: 
since we assumed that k and \h\ are statistically indepen- 
dent, A.K c {\h\) is determined solely by the variance a\ 
of the threshold distribution around the mean threshold 
I h I = \h\ - the smaller this variance is, the more peaked 
this distribution is around \h\ — \h\, and the less it hence 
differs from the homogeneous distribution. Since we as- 
sumed that (in the inhomogeneous case) thresholds are 
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Poisson distributed around \h\, we directly conclude 

AK c (\h\) ~ a 2 h = \h\. (25) 

For arbitrary threshold distributions that are statistically 
independent from the networks' degree distribution with 
variance a\ , we make the ansatz 

<7it=4 + <rl (26) 

for the total variance af ot . Using the same Gaussian ap- 
proximation as above for the homogeneous case, one can 
show that 



K c {\h\) 



21n|/i| 



(27) 




for networks with inhomogeneous thresholds distributed 
around an average absolute threshold \h\ (for details, cf. 
appendix C). This implies that, in the limit \h\ — > oo, all 
networks with Poissonian distributed connectivity and 
threshold distributions with a variance which obeys the 
scaling relation: a\ ~ \h\@ with < j3 < 2, follow the 
universal asymptotic scaling relation 

Jr 



K c {\h\) = 



2\xy\hV 



(28) 



as it is shown in appendix C. This means that in all these 
cases, the asymptotic scaling for \h\ — * oo is dominated 
by by the scaling behavior of the maximum of the damage 
propagation function p s (k, |/i|), with an exponent a = 2. 
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FIG. 9: /i'c ( | ft. | ) for homogeneous thresholds (+) and Poisson 
distributed thresholds with the same average \ h\ (X), annealed 
approximation. The solid line is the asymptotic scaling ob- 
tained from Eq. (|15jl . For inhomogeneous \h\, the critical line 
is systematically below K c of networks with homogeneous \h\. 
Inset: The difference | Aii'cd ft| ) | between both curves grows 
only linearly in \h\, confirming that the asymptotic scaling in 
the limit \h\ — > oo, is the same in both cases. 



FIG. 10: K c (/3, \h\) for networks with threshold distribu- 
tions following discretized Gaussian distributions with differ- 
ent variances Var(\h\) — \h\^ (for details, see text). One 
clearly appreciates that the larger the variance of the thresh- 
old distribution, the more the curves K c ((3, \ h\) are below the 
critical line of networks with homogeneous thresholds (blue 
solid line); in the limiting case /3 = 1.95 w a (yellow trian- 
gles), K c scales almost linearly with \h\. Inset: differences 
| AK c (/3, \h\)\ to the critical line of RTN with homogeneous 

thresholds scale ~ \h\^° with f3 < f3 e < a (power law fits and 
dashed line with slope a shown for comparsion); this implies 
asymptotic convergence to the universal scaling function Eq. 
(1281) in the limit \h\ —> oo for all cases shown here. 



p 




0.5 


0.533 ±0.009 


1.0 


1.099 ±0.004 


1.2 


1.327 ±0.004 


1.5 


1.732 ±0.004 


1.8 


1.942 ± 0.003 


1.95 


1.975 ±0.004 



TABLE I: Scaling exponents (3 e , as obtained from fits of 
AK C ~ jftj' 3 ", as a function of /3. 



Let us now confirm this finding for a different class of 
threshold distributions. Since in a Poissonian the vari- 
ance is not a free parameter, we now instead choose a 
discretized Gaussian distribution, i.e. 



P(\h\) = 



Z 



with 



Z 



and variance 



1 



\h\=0 



(JhV^TT 



e 2 



i(N-|fc|) 2 K 



(29) 



(30) 



a 2 h = \hf, 0£[O,a). 



(31) 
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The factor Z ensures that the probabilities are normal- 
ized in the interval [0, |/i| m ], where \h\ m denotes the cut- 
off of the threshold distribution. Figure fTOl compares the 
scaling functions i£" c (|/i|) for different values of (3 to the 
asymptotic case of homogeneous networks. Obviously, 
for finite \h\, increased variance of the threshold distri- 
bution substantially lowers the critical connectivity; in 
the limiting case (3 ~ a, K c grows only linearly with \h\. 
For (3 < a, we find that the deviation from the scaling 
behavior of RTN with homogeneous thresholds scales as 



10000 



AK C oc \hf 



(32) 



Table 1 compares (3 and j3 e (as obtained from fits of AK C ; 
in all cases, we have [3 e > (3, which is a discretization 
effect, but still (3 e < a. Hence, it follows that 



hm K ^£ll ]hl) 



= lim 

\h\ — >oo 

= 1 — const. 
= 1 



K£(\h\) - AK c (P,\h\) 



KH\h\) 

lim \hf'- a (33) 

\h\ — >oo 



for (3 e < a, i.e. in this case all scaling functions K c ((3, \h\) 
for \h\ — » oo indeed asymptotically converge to the same 
universal scaling function, as given by Eq. (|28[) . 

Let us now have a closer look at the scaling behav- 
ior of the intersection points .fQ(|/i|), as introduced in 
the last paragraph of subsection B. Let d h (K, \h\) be the 
expected damage in networks with homogeneous thresh- 
old, and d l (K, \h\) the expected damage in networks with 
inhomogeneous thresholds; then 

d h {K d {\h\),\h\)-d\K d {\h\),\h\) = (34) 

\h\ = \h\ (35) 

are the defining equations for Kd(\h\). Notice that for 
K < K dl the randomness introduced by inhomogeneous 
thresholds actually increases the probability for damage 
spreading, whereas for K > Kd, it is decreased. Equation 
(|34|) . under condition Eq. (|35|) . can be solved numerically 
for not to large \h\. Further, one can derive the asymp- 
totic scaling in the thermodynamic limit by application 
of the Gaussian approximation for the damage propaga- 
tion function (for details, cf. appendix D), showing that 



lim K d (\h\) = h 2 - \h\. 



\h\- 



(36) 



Fig. [TT] demonstrates that Kd{\h\) approaches this 
asymptotic scaling already for considerably small \h\, in- 
dicating that Kd(\h\) is characterized by the same uni- 
versal scaling exponent a = 2 as K c (\h\). Notice, how- 
ever, that the asymptotic scaling law for Kd obeys a 
purely algebraic relation, whereas K c has a dependence 
~/i 2 /ln|/i| (Eq. El- 
Let us briefly compare the scaling behavior of RTN 
with non-zero thresholds, as discussed above, to Ran- 
dom Boolean Networks (RBN). Obviously, increasing \h\ 




FIG. 11: Scaling behavior of Kd(|/l|) as a function of \h\, dou- 
ble logarithmic plot. The dashed line highlights the asymp- 
totic scaling (Eq. (|36p h 



biases the output states of network nodes (for the sys- 
tems discussed in this paper, it increases the probability 
to have an output state <x; = —1). Biased RBN obey the 
scaling relationship 0] 



Kr 



1 



2p(l-p) 



(37) 



To compare this relationship to the asymptotic scaling 
for RTN in the limit of large \h\, we have to consider 
the limit p — > 1. One can show that, in this limit, the 
scaling function Eq. ([37|) logarithmically approaches the 
asymptotic scaling 



p- 



21np 



(38) 



This shows that \h\ plays the same role as the bias pa- 
rameter p in RBN, and that both classes obey the same 
scaling in the limit p — > 1 and \h\ — > oo, respectively. 
However, there are also substantial differences between 
both classes of systems, that come into play when \h\ is 
small (when p is close to 1/2). In particular, while RBN 
in this limit still obey the simple scaling relationship Eq. 
(|37p). the critical connectivity K c of RTN is derived from 
the complex dependence of Eq. (fl~4"|) . This difference is 
due to the fact that, in RTN, local damage propagation 
strongly depends on the in-degree of nodes (cf. Eq. [6] and 
IH]), while it is independent from the in-degree in RBN for 
k > 0. In the limit of sparsely connected networks (i.e. 
small \h\ and K c ), this leads to much stronger finite size 
effects in RTN than in RBN. Furthermore, in this limit 
also the absolute values of K c in RTN are considerably 
below those of RBN H Hj| . 

Finally, let us remark on the existence of the character- 
istic connectivity K d . As shown above, Kd is defined for 
RTN with arbitrary \h\, in particular, it exists in the limit 
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\h\ — > oo, with a well-defined asymptotic scaling. For bi- 
ased RBN, the corresponding limit is given by p — > 1 (or, 
equivalently, p — > 0). Obviously we can in principle as- 
sign variable (inhomogeneous) biases pt to different RBN 
nodes such that the average bias is equal to p. How- 
ever, because p is a probability and hence < p < 1, 
the variance has to vanish in the limit p — > 1 (p — > 0) 
to yield a proper average bias. Since Kd is defined by 
comparing networks with diverging variance of the order 
parameter \h\ (or p, respectively) with the corresponding 
networks with vanishing variance and the same average 
| ft | (or p, respectively), this implies that Kd is not de- 
fined for RBN in the limit of large bias p — ► 1, which 
corresponds to |ft| — * oo in RTN. Hence, Kd constitutes 
a truly novel, not previously known concept, yielding a 
new characteristic connectivity which is well-defined only 
for RTN. 

It is interesting to notice that the dependence of K c , 
as well as of Kd on |ft| is clearly super-linear even for 
considerably small |ft|; this has profound consequences 
for algorithms that evolve RTN towards (self-organized) 
criticality by local adaptations of both thresholds and the 
number of inputs a node receives from other nodes [33 | . 
In particular, it can be shown that co-evolution of net- 
work dynamics and thresholds/in-degrees leads to strong 
correlations between |ft| and k. To approach this type 
of problem analytically, we will now extend our analysis 
in this direction, first, In the next section, we will show 
that even weak correlations between k and |ft| can lead to 
a transition from sub-critical to super-critical dynamics 
(and vice versa), while keeping the average connectivity 
K and the average absolute threshold |ft| constant. 

Generation of correlations: 
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FIG. 13: Combined density p(ki„, \h\) for three different pos- 
itive values of the correlation parameter c, from top to bot- 
tom: c — 0, c — 0.2 and c = 0.9. Dark gray indicates a high 
probability density. The diagonal structure of p(kin, \ h\) for 
c = 0.9 (lower panel) indicates emergence of strong positive 
correlations between ki n and \h\. 



FIG. 12: Schematic illustration of the algorithm applied to 
generate local (anti-)correlations between in-degree ki n and 
(absolute) threshold \h\. Arrows symbolize inputs from other 
nodes, boxes symbolize node thresholds (one box corresponds 
to \h\ — 1, two boxes to \h\ = 2, and so on). For details of 
the algorithm, please refer to the text. 



D. Effect of correlations between k and h 

So far, we assumed that node degree and node thresh- 
olds are totally uncorrelated; while this matches well the 
" maximum disorder" assumption used in random ensem- 
ble based approaches as, e.g., the annealed approxima- 
tion, this might be a quite unrealistic constraint for many 
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FIG. 15: Average damage d(c) as a function of |c|, for corre- 
lated ki n and (+) and anti-correlated ki n and (X), with 
K — 6.15 for c > networks, A' = 5.8 for c < networks and 
\h\ = 2.5 in both cases. Numerical data where obtained from 
ensemble averages over Z — 5 • 10° randomly generated RTN 
with N — 1024 nodes for each data point. Solid curves are 
the corresponding results of the corrected annealed approxi- 
mation, while the dashed-dotted curve shows the uncorrected 
result for correlated i n k and \h\, and the dashed curve the un- 
corrected result for anti-correlated k in and \h\, respectively. 
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\h\ 

FIG. 14: Combined density p(ki n , \h\) for three different neg- 
ative values of the correlation parameter c, from top to bot- 
tom: c = 0, c — —0.1 and c = —0.9. Dark gray indicates a 
high probability density. The inverted diagonal structure of 
p(ki n , \h\) for c = —0.9 (lower panel) indicates emergence of 
strong anti-correlations between fc;„ and \h\. 

real world networks. Indeed, one can show that even in 
a simple evolutionary algorithm that couples both the 
adaptation of node thresholds hi and in-degree ki to a 
local dynamical order parameter, strong correlations be- 
tween both quantities emerge spontaneously (32J. Hence, 
it is an interesting question to ask whether correlations 
(or anti-correlations) between h and k may induce a tran- 
sition from sub-critical to super-critical networks (or vice 



versa), while we keep K and \h\, and network topologies 
constant. 

Let us first formulate an algorithm that generates 
correlations (anti-correlations) between k and \h\. For 
this purpose, a parameter c G [0,1] is introduced which 
parameterizes the probability that k and \h\ are locally 
correlated (anti-correlated). The topology-generating 
algorithm then reads as follows (compare also Fig. [T2")) : 

1) Generate a random, directed network with Pois- 
son distributed k and Poisson distributed \h\ with 
average connectivity K and \h\ for all sites. 

2) Select a pair of sites i < N and j < N at random, 
c > 0: exchange the sites' thresholds if fcj n (i) > kin(j) 
and \h\(i) < \h\(j), or vice versa, c < 0: exchange the 
sites' thresholds if ki n (i) < ki n (j) and \h\(i) < \h\(j), or 
vice versa. 

3) Go back to 2) and repeat the algorithm for c x P rnax 
steps, where P m ax is a pre-defined maximum number of 
correlated pairs. 

Obviously, increasing the parameter c € [0, 1] increases 
correlations (anti-correlations) between ki n and h. If wc 
repeat this algorithm Z times for fixed c, we can gen- 
erate a random ensemble of Z correlated/anti-correlated 
networks, and investigate damage spreading on these net- 
works. The ensemble-averaged probability p(fcj n , \h\) to 
have a site with fej„ inputs and threshold = \h\ then 
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FIG. 16: Difference Ap(k, \ h\) between real and effective com- 
bined densities (see text), for three different values of positive 
c. Dark gray indicates a negative deviation, light gray positive 
deviation, a medium intensity refers to zero deviation. 



is defined as 



p{hn, \h\) 



T,j=l n j( k in, \ h \) 



Z-N 



(39) 



where rij{ki ni \h\) is the number of sites with ki n inputs 
and threshold \hi\ = \h\ in the jth random network. Fig- 
ure [13] demonstrates the correlating effect of the algo- 
rithm on the average probabilities p(ki n , \h\) for ensem- 
bles of 10 5 randomly generated networks, for the case 
c > 0, with P max = 10 4 . For c = 0, clearly no cor- 




- ft 



K 



FIG. 17: Schematic illustration of the naive mean-field as- 
sumption implicit to the annealed approximation: Choosing 
a random ensemble of K nodes and inverting one input of 
each of them (left panel, black circles refer to inverted states) 
on average yields the same damage as perturbing a randomly 
chosen site with K outputs and investigating the resulting 
damage at the output nodes (right panel); in both panels, X 
means "damaged states". This assumption is violated for cor- 
related networks, even if the generating algorithm explicitly 
correlates only in-degree and thresholds. 



relations are present, and the combined density simply 
represents the independent superposition of the two un- 
derlying Poisson distributions. With increasing c, corre- 
lations gradually emerge, and for c = 0.9 the resulting 
distribution clearly exhibits a diagonal structure. Figure 
[T4l demonstrates the corresponding effect for c < 0, i.e. 
anti-correlated topologies. 

Let us now investigate how these correlations affect 
damage propagation. To apply the annealed approxima- 
tion, we now have to calculate the average probability 
for damage propagation (in a finite network of size N) 
according to 

\h\ m N 

<p.)(K,\h\,c)= X) P(kin,\h\)p s (k + l,\h\), (40) 

\h\=0 k=\h\ 



with the normalization conditions 

\h\ m N 

£ Yl P{hnM) = l 

|fc|=0 k=\h\ 



(41) 



and 



\h\ m N 

E \Hp{kin,\h\) 

\h\=0 k=\h\ 



(42) 



where \h\ m is the maximal absolute threshold observed 
(cutoff); correlations enter via the probabilities p(ki n , \h\) 
to observe a node with in-degree fcj„ and absolute thresh- 
old \h\. 

Figure [15] compares the numerically observed damage 
d (for ensembles of randomly generated networks) one 
time step after a_ one-bit perturbation to the expected 
damage (p s )(K, \h\,c) ■ K, as predicted by the annealed 
approximation (dashed curves). In both cases, for corre- 
lations and anti-correlations, the so-obtained theoretical 
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curves are obviously wrong both with regard to quanti- 
tative matching of the numerical data, as with regard to 
the predicted trend: while in the numerical experiment 
a decrease of d with increasing c > is observed, i.e. a 
transition from super-critical (chaotic) to sub-critical (or- 
dered) dynamics, the annealed curve predicts a strong in- 
crease. Corresponding observations (with opposite signs) 
are made for the case of anti-correlations. We conclude 
that there must be an additional effect present which is 
not captured in our naive mean-field model. To iden- 
tify the origin of this deviation, one has to compare the 
distribution p(ki n , \h\) which is observed for the outputs 
of pertubed sites with the original distribution p(ki n , \h\), 
which is averaged over the whole topology. Figure 1161 
shows the deviation Ap(foj n , \h\) := p(ki n , \h\)—p(ki n , \h\) 
between both distributions. While for small ki n and \h\ 
negative deviations are found, for larger ki n and \h\ de- 
viations have a positive sign. One can easily understand 
the source of this effective bias, when one thinks of the 
correlation-generating mechanism (Fig. 112ft : if we pick, 
by chance, within a correlated network a site i with small 
ki ni it will probably also have a small \h\, and its outputs 
will probably have larger ki n and \h\ than site i. Since 
both ki n and \h\ are bounded from below, this leads to 
a systematic bias to observe larger ki n and \h\, at the 
expense of smaller, for the outputs of perturbed nodes. 
The opposite effect obviously holds for sites with large 
km and \h\. While for c = this effect is very small, 
it becomes dominant for c — > 1 due to the strong asym- 
metry of the diagonal distribution. If we correct the an- 
nealed approximation to include this bias, by replacing 
p(kin, \h\) with p(kin, \h\) in Eq. (|40|) . we find that the 
resulting corrected annealed curves much better match 
the numerically observed damage (still, there are slight 
discrepancies for large c, which are due to finite ensemble 
sizes) . 

The result of this study demonstrates that the an- 
nealed approximation has to be used with extreme care, 
when topological correlations are present. Although the 
applied algorithm explicitly correlates only in-degree and 
thresholds, consideration of these correlations only by us- 
ing a combined density p{k% n , \h\) for averages over topol- 
ogy and dynamics leads to wrong predictions. In addi- 
tion, one has to consider systematic bias effects between 
perturbed sites and their outputs, which arise as a side 
effect of the correlating algorithm. This shows that the 
naive mean-field assumption inherent to the annealed ap- 
proximation, as demonstrated in Fig. 1171 is violated even 
in the simplest case of correlations between in-degree and 
thresholds. Instead, complete information about topol- 
ogy, including the output side, has to enter statistics, and 
the order of averages becomes important. 



IV. DISCUSSION 

An increasing number of studies is concerned with the 
propagation dynamics of perturbations and/or informa- 



tion in complex dynamical networks. Discrete dynam- 
ical networks, in particular Random Boolean Networks 
(RBN) and Random Threshold Networks (RTN), consti- 
tute an ideal testbed for this type of question, since they 
are easily accessible for both computational methods and 
the tool boxes of statistics and combinatorics. Often, it 
is found that damage/information propagation strongly 
depends on the type of inhomogeneities present in net- 
work wiring. Several studies focus, for example, on the 
effect of scale-free degree distributions [1^, [26| . Typically, 
these studies employ mean-field methods and hence rep- 
resent, in a sense, strongly idealized models, since they 
derive results that strictly hold in the thermodynamic 
limit only. 

Consequently, a second line of research concentrates 
on modification of damage propagation due to finite-size 
effects, which play a decisive role in many real- world net- 
works. Recently, it was shown that weakly perturbed, fi- 
nite size RBN and RTN show pronounced deviations from 
the annealed approximation [il| . Fronczak and Fron- 
czak showed that these deviations can be explained by 
inhomogeneities and emergent correlations found at the 
percolation transition [33j . however, their study is cur- 
rently limited to undirected networks. In this context, 
the system discussed in our paper constitutes a comple- 
mentary approach: it allows to introduce dynamical in- 
homogeneity of network units, without otherwise altering 
network topology. While this type of dynamical diversity 
certainly plays an important role in many real- world net- 
works, it is neglected by most researchers. Let us now 
briefly summarize the main results of our study. 

We studied damage propagation in Random Thresh- 
old Networks (RTN) with homogeneous and inhomoge- 
neous negative thresholds, both analytically (using an 
annealed approximation) and in numerical simulations. 
We derived the probability p s (k, \h\) of damage propa- 
gation for arbitrary in-degree k and (absolute) threshold 
\h\ (Eqs. (5)-(8)), and, from this, the corresponding an- 
nealed probabilities (p s ) (Eq. (10) and Eq. 12)) and the 
expected damage d (Eq. (11)), for both the cases of ho- 
mogeneous and inhomogeneously distributed thresholds. 
On these grounds, we investigated the scaling behavior 
of the critical connectivity K c as a function of \h\. Using 
a mean field approximation, a simplified scaling equa- 
tion for the logarithm of the average damage was derived 
(Eq. (14)), and applied to derive the critical line lf c (|/i|) 
(Fig. 6). It was shown that this function exhibits a super- 
linear increase with \h\, which asymptotically approaches 
a unique scaling law if c (|/i|) ~ h 2 /(2 In \h\) for large \h\ 
(Eq. (18) and Fig. 7). However, convergence against this 
asymptotic scaling is very slow (logarithmic in \h\), which 
indicates that finite size effects are very dominant, and 
cannot be neglected for realistically sized networks. We 
presented evidence that this asymptotic scaling is uni- 
versal for RTN with Poissonian distributed connectivity 
and threshold distributions with a variance that grows 
slower than h 2 , for both the cases of Poisson distributed 
thresholds (Fig. 8) and thresholds distributed according 
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to a discretized Gaussian (Fig. 9). Interestingly, inho- 
mogencity in thresholds, meaning that each site has an 
individual threshold \hi\ drawn, e.g., from a Poisson dis- 
tribution with mean \h\, increases damage for small aver- 
age connectivity K, when compared to homogeneous net- 
works with the same average threshold \h\ — h, whereas 
for larger K with K > Kd, damage is reduced. This es- 
tablishes a new characteristic connectivity K"<j(|/i|) with 
Kd > K c , that describes the ambivalent effect of thresh- 
old inhomogeneity on RTN dynamics. We showed that 
Kd (\h\) asymptotically converges against a unique scaling 
law Kd ~ h 2 in the limit \h\ — > oo. The scaling of Kd was 
compared to the corresponding case of random Boolean 
networks (RBN) with inhomogeneously distributed bias, 
parameterized in terms of a bias parameter 1/2 < p < 1. 
It was shown that Kd is not defined for RBN in the limit 
p — ► 1, which corresponds to \h\ — > oo in RTN. Hence, 
Kd constitutes a truly novel, not previously known con- 
cept, yielding a new characteristic connectivity which is 
well-defined only for RTN. 

Last, we introduced local correlations between in- 
degree ki n of network nodes and their (absolute) thresh- 
old | ft, |, while keeping all other network parameters con- 
stant. We found that even small positive correlations can 
induce a transition from supercritical (chaotic) to sub- 
critical (ordered) dynamics, while anti-correlations have 
the opposite effect. It was shown that the naive mean- 
field assumption typical for the annealed approximation 
leads to false predictions in this case. Even in the sim- 
plest case, where only in-degree and (absolute) threshold 
are correlated, complete information about topology, in- 
cluding the output side, has to enter statistics, and the 



order of averages becomes important. 

To summarize, dynamics of damage (or information) 
propagation in RTN with inhomogeneous thresholds and 
Poisson distributed connectivity shows both similarities 
and differences, when compared to networks with ho- 
mogeneous thresholds: similarities manifest themselves 
in common universal scaling functions for both K c and 
Kd, whereas differences show up in the opposite effects 
of threshold inhomogeneity for small and large K. Dif- 
ferences become even more prominent in networks that 
are characterized by correlations between in-degree and 
thresholds. In this case, the annealed approximation has 
to be used with extreme care. Many dynamical systems 
in nature, that can be described as complex networks, 
exhibit considerable variation of activation thresholds 
among the elements they consist of, however, these varia- 
tions are often neglected (e.g., in Boolean network based 
models of gene regulation networks). Our results indi- 
cate that, while general characteristics as, for example, 
the scaling behavior of critical points, may be conserved 
in approxmations of this type, inhomogeneous thresholds 
can strongly impact the details of network dynamics, and 
hence should be taken into account in models that aim to 
give a realistic description of the dynamics of, e.g., gene 
regulation networks. 
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[34] We restrict ourselves to negative (or zero) thresholds, to 
ensure that the 'default state' of a network site i, i.e. 
when its inputs sum to zero, is to be 'inactive' [en — —1), 
which naturally excludes positive thresholds. 
[35] Other authors define sgn(0) = +1, however, for sym- 
metry reasons update dynamics is not affected by either 
choice. If we interpret the state <7; = — 1 as 'inactive' and, 
correspondingly, +1 as 'active', our choice appears to be 
more natural: the default state of a network site is to be 
'inactive', unless it receives activating inputs from other 
sites. 

[36] To obtain accurate results, one has to consider networks 
sizes N K, and adjust the upper limit of the sum in 
(|10[) accordingly. Since a small step size AK has to be 
applied iteratedly to identify K c , this becomes computa- 
tionally very costly. 



Hence, the damage propgation probability follows as 



,(k,\h\) = k- 1 -2-( fc+1 ) • 



(k + \h\ + l)- f fc+ |h|+i ) 



-(k-\h\ + l)- fc _ N 



l+i 



2-(*-i)(fc- 1)! 



[{k + \h\-l)/2]l[(k-\h\-l)/2}\ 



(A5) 
(A6) 



,-(fc-i) ( ( fc 1 



= 2 ~ (k ~ l \^k)- (A7) 

2) k — \h\ even: Here, we have as necessary conditions 
k+ -k--\h\ = Q (A8) 



\h\ 



(A9) 



In the first case, only the reversal of negative spins is 
effective, whereas in the latter case the same holds for 
positive spins. We have 



(A10) 



APPENDIX A: DERIVATION OF p s (k,\h\) 

In this section, we provide a derivation of the local 
damage propagation probability p s (k, \h\. 

Consider a network site i with k inputs; k + of these 
have positive sign, fc_ negative sign, hence, k + + fc_ = k. 
We no derive the conditions under which a inversion of 
one input spin at time t leads to a switch of the output 
of site i at time t + 1. 

1) k — \h\ odd: From Eqs. [T] and Git is easy to see 
that input-spin flips produce " damage" only if one of the 
following conditions holds: 



in the first case and 



\h\ = l 



k+ - fc_ - \h\ = -1. 



(Al) 



(A2) 



In case ED only the reversal of positive spins is effective, 
whereas in case IA21 only the reversal of negative spins 
has an effect. We have 



k+ = 



in the first case and 



k+\h\ + l 



k -\h\ + l 



(A3) 



(A4) 



in the second case. There is a total number of k ■ 2 



possible spin configurations, of which 



fill condition IA3I and ( 



.(k+\h\+l)/2) 



ful- 



(k-\h\ + l)/2) 



fulfill condition IA4 



fc+ \h\ + 2 



(AH) 



in the second case. There is a total number of k ■ 2 k pos- 
sible spin configurations, of which {^ k _^y 2 ) fulfill con- 
dition [XTU] and ((fc + |/ l | f + 2)/2) fulfill condition I Al 11 Hence, 
the damage propgation probability follows as 



Ps(k,\h\) = fc" 1 ^-^ 1 ) 



(*-W)-(jHh|) 



+ (k+\h\+2)- (k+jhj+2 
V 2 

2 -(fc-i)(fc- l)! 
[(k-\h\-2)/2}\[(k + \h\)/2]\ 

- z i k+\h\ 



(A12) 
(A13) 
(A14) 



APPENDIX B: DERIVATION OF THE SCALING 
EQUATION 

For RTN with Poisson distributed in- and out-degree, 
the critical line is given by the condition 

d(t + 1) = ( Ps )(K c (\h\), \h\) ■ K c (\h\) = 1. (Bl) 



with 



N 



(p s )(K, \h\) = e~ R J2 ^-Ps(k + 1, \h\). (B2) 



k=\h\ 



k\ 
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Instead of averaging over the ensemble of all possible net- 
work topologies as in Eq. (|B2|) , we now make an explicit 
mean field approximation, and consider a "typical" net- 
work node with k rs K inputs. Consequently, we approx- 
imate 



(p s )(K,\h\)^p s ([K\,\h\), 



(B3) 



where [-J denotes the floor function. In the limit of large 
K and \h\, the difference between the damage propaga- 
tion probabilities for even and odd k vanishes, i.e. we 
can set 



<P.>(^|fc|)«2-(L*J-«f (L J |J i: 



and hence 



d(K,\h\) = K-2-^i 



[K\ 

2 



(B4) 



(B5) 



without loss of generality. 

Using the Stirling approximation n\ « n n e~ n \phrn, 
dropping the floor function (since we now consider a func- 
tion of real- valued variables only) and taking logarithms, 
we obtain 

In [d(K, \h\)\ w InK - bx2- K + Z x - Z 2 - Z 3 (B6) 
with 



In 



Z x = \n[K"< 
K-\h\ 



K„-K 



V2ttK], 



n(K-\h\) 



and 



Zz = In 



K + \h\ 



n(K+\h\) 



Summing out the logarithms in Z%, Z2 and Z3, one 
realizes that all terms linear in K drop out, resulting in 



In [d(K, 



InK 
K-\h\ + l 



K - - ) In A 



2 

K+\h\ 



1 



\n{K-\h\) 

In [K+\h\) + C{B7) 

with C = In ^ 2 /-7rj . Using some simple algebra and 
approximating \h\ + 1 w \h\, this can be reformulated as 

\n[d(K,\h\)} « lnK--{ln(K 



-Kin 
h\h\ln 



(K+\h\)(K-\h\) 



K 2 



K+\h\ 
K - \h\ 



C. 



(B8) 



This leads to the final result 
1 



]n[d(K,\h\)] 



2 

- \h\ln 



\nK -K -\n 
K+\h 



[K-\h\ 



C. 



(B9) 



APPENDIX C: ASYMPTOTIC SCALING OF K c 

Let us now derive the asymptotic scaling behavior of 
the critical connectivity K c (\h\). We start with the case 
of homogeneous thresholds, and then generalize to inho- 
mogeneous thresholds. 

First, we note that the right hand-side of Eq. (|B5|) has 
the form of a Binomial distribution 



P{n,k) 



n jn — h 



p q 



(CI) 



with p = q = 1/2, n = [K\ and k = ([K_\ + \h\)/2, 
multiplied with a prefactor K . In the limit K — > 00 and 
I ft. I — > 00, we can replace the Binomial distribution with 
a Gaussian and drop the floor function, i.e. 



d(K, \h\) = K-C n exp 



K+\h\ 
2 



This simplifies to 
d(K, \h\) 



= K- 



2ttK 



exp 



2K 



(C2) 



(C3) 



with the normalization constant C n = ^J\j (2wK) and 
variance a 2 = K. 

In the case of inhomogeneous thresholds, we can still 
use this approximation, however, the variance of the 
threshold distribution adds to the variance of the dam- 
age propagation function of the homogeneous case. This 
implies that we have to replace K with K+a^, and hence 



d(K,\h\) = 



K + o-l 



2n(K + a 2 h ) 



exp 



2(K + a 2 ) 



(C4) 



To obtain the criticality condition, we take logarithms 
and set the result to zero, leading to 



In [K c + 0-1] -\ In [2n(K c + a 2 h )] - —A- 

2 2{R c + a h ) 



This simplifies to 



{K c 



a?) In 



K r . 



2tt 



= 0. 

(C5) 

(C6) 



To solve this equation with respect to K Cl we make the 
ansatz 

h 2 



K c + at 



21n|/il 



(C7) 
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Inserting for K c + a\ into Eq. (|C6|) , we obtain 



■In 



21n|/i| 



47rln h\ 



= h 2 il- 



ln [47T In \h] 
2hi\h\ 



(C8) 
(C9) 



Since the second term in the bracket vanishes logarith- 
mically for | h | — > oo, we have verified that Eq. (|C7|) 
yields the correct asymptotic scaling. Consequently, the 
asymptotic scaling of the critical line for large \h\ is given 
by 



2 In /i 



'ft- 



(CIO) 



However, notice that the convergence is very slow, as can 
be appreciated from the logarithmic finite-size term in 
Eq. (|C9p . In particular, we conclude that the asymptotic 
scaling for networks with homogeneous thresholds, i.e. 
\h\ = const, and er/j = is given by 



21n ' 



(Cll) 



Let us now prove that this scaling is universal for \h\ — > 
oo for all threshold distributions possessing a variance 
ai ~ \h\ a with < a < 2. In this case, we have 



1- Pj, 

Hhoo AT£ om (|/l|) 



lim 



K c hom (l/ll) 



i^ om (l/ll) 



— (C12) 



= 1 — lim 



|fc|-»ao K% om {\h\) 



(C13) 



= 1 — lim 



21n|/i| 



|/i|->oo 



2-a 



(C15) 



Since we assumed < a < 2, the limit in Eq. (|C15[) van- 
ishes, and hence the asymptotic scaling equation IC11I is 
indeed universal for this class of threshold distributions. 



APPENDIX D: ASYMPTOTIC SCALING OF K d 



networks. Let us further assume that thresholds are Pois- 
sonian distributed, i.e. a\ — \h\. If we apply the same 
Gaussian approximation as in section[Cj these conditions 
lead to 



-h 2 /(2K d ) 



-h 2 /(2(K d + \h\)) 



V^iQ ^{K d + \h\) 

Taking logarithms and reordering, this reduces to 



In 



K d +\h\ 
K d 



K d K d + \h\ 







(D3) 



(D4) 



Linearization of the first term leads to the approximation 



K~ d ~ K~ d + K d + \h\ 



(D5) 



Solving this equation for K d finally yields the asymptotic 
scaling 



K d (\h\) « h 2 - \h\, 
i.e. K d scales quadratically with \h\. 



(D6) 



APPENDIX E: POWER-LAW APPROXIMATION 

OF K c (\h\) FOR FINITE |fe| 
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The characteristic connectivity K d is defined by the 
conditions: 



\h\ = \h\, 



(Dl) 



where \h\ is the (constant) threshold of a homogeneous 
network, and \h\ the average threshold of a corresponding 
network with inhomogeneous thresholds, and 

d h (K d (\h\), \h\) - d\K d {\h\), \h\) = 0, (D2) 

where d h is the expected damage for homogeneous net- 
works, and d % is the expected damage for inhomogeneous 



FIG. 18: Solutions of Eq. (JE3J for (from the left to the right) 
a — 1.6, q = 1.7, a — 1.8 and a — 1.9. Projections of the 
maximum on the |/i|-axis (as indicated by arrows) yield the 
corresponding values of \h\ c at which the approximations are 
optimal. 

In this section, we first describe how to identify nu- 
merically candidate solutions (power-laws) 



K c {\h\)*a(\h\)'\h\ aW - 



(El) 



that optimally approximate Eq. for finite (critical) 
\h\c 
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We start with a fixed a £ [1.6, 2) and define 



F(y) 



(|/i|+l)ln 









■f In y — y ■ In 


1 - 









lV-\h\ 



C 



(E2) 



with y — a ■ \h\ a . One can show that, for any finite a 
and a, F(y) has a maximum at a finite value |/i| maa :- We 
know that K c is a monotonically increasing function of 
\h\, and intend to optimize the power-law approximation 
exactly at K c . Hence, we have to vary a such that 



maxF(?/)| Q = 0. 



(E3) 



Projection of the maximum on the |ft,|-axis then yields 
the corresponding threshold values | c (ct) at which the 
approximation for the given a is optimal (Fig. [T5]). In- 
version of this relation allows us to plot the corresponding 
values of the function a(|/i|) (Fig. [5]). 

Last, let us estimate the asymptotic scaling of a(|/i|). 
If we apply the asymptotic scaling relation for K c derived 



in section [Cj we can approximate 
h 2 



a(\h\) ■ \h\ a ^l 



(E4) 



2\n\h\ 

Taking logarithms, this yields 

2 ln\h\ - In 2 - In In \h\ = lna(\h\) + a(\h\) hi \h\. (E5) 

We now consider variations of a only, i.e. we fix a with 
respect to \h\. Taking the derivative with respect to \h\ 
on both sides of the equation and solving for a then yields 



a(\h\) 



1 



ln|ft| 



(E6) 



Inserting this result into Eq. (|E5p . we finally obtain the 
estimate 



a(\h\) 



21n|/i| 

for the proportionality constant a. 



(E7) 



